from eofs.xarray import Eof
import xarray as xr
import numpy as np
import datetime as dt
import matplotlib.pyplot as plt #用来绘图的包
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.ticker import (LongitudeFormatter, LatitudeFormatter, LatitudeLocator, LongitudeLocator)
import cmaps
file = xr.open_dataset('./mydata/sst_month_mean.nc')
file
<xarray.Dataset>
Dimensions: (lat: 89, lon: 180, time: 2019)
Coordinates:
* lat (lat) float32 88.0 86.0 84.0 82.0 80.0 ... -82.0 -84.0 -86.0 -88.0
* lon (lon) float32 0.0 2.0 4.0 6.0 8.0 ... 350.0 352.0 354.0 356.0 358.0
* time (time) datetime64[ns] 1854-01-01 1854-02-01 ... 2022-03-01
Data variables:
sst (time, lat, lon) float32 ...array([ 88., 86., 84., 82., 80., 78., 76., 74., 72., 70., 68., 66.,
64., 62., 60., 58., 56., 54., 52., 50., 48., 46., 44., 42.,
40., 38., 36., 34., 32., 30., 28., 26., 24., 22., 20., 18.,
16., 14., 12., 10., 8., 6., 4., 2., 0., -2., -4., -6.,
-8., -10., -12., -14., -16., -18., -20., -22., -24., -26., -28., -30.,
-32., -34., -36., -38., -40., -42., -44., -46., -48., -50., -52., -54.,
-56., -58., -60., -62., -64., -66., -68., -70., -72., -74., -76., -78.,
-80., -82., -84., -86., -88.], dtype=float32)array([ 0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20., 22.,
24., 26., 28., 30., 32., 34., 36., 38., 40., 42., 44., 46.,
48., 50., 52., 54., 56., 58., 60., 62., 64., 66., 68., 70.,
72., 74., 76., 78., 80., 82., 84., 86., 88., 90., 92., 94.,
96., 98., 100., 102., 104., 106., 108., 110., 112., 114., 116., 118.,
120., 122., 124., 126., 128., 130., 132., 134., 136., 138., 140., 142.,
144., 146., 148., 150., 152., 154., 156., 158., 160., 162., 164., 166.,
168., 170., 172., 174., 176., 178., 180., 182., 184., 186., 188., 190.,
192., 194., 196., 198., 200., 202., 204., 206., 208., 210., 212., 214.,
216., 218., 220., 222., 224., 226., 228., 230., 232., 234., 236., 238.,
240., 242., 244., 246., 248., 250., 252., 254., 256., 258., 260., 262.,
264., 266., 268., 270., 272., 274., 276., 278., 280., 282., 284., 286.,
288., 290., 292., 294., 296., 298., 300., 302., 304., 306., 308., 310.,
312., 314., 316., 318., 320., 322., 324., 326., 328., 330., 332., 334.,
336., 338., 340., 342., 344., 346., 348., 350., 352., 354., 356., 358.],
dtype=float32)array(['1854-01-01T00:00:00.000000000', '1854-02-01T00:00:00.000000000',
'1854-03-01T00:00:00.000000000', ..., '2022-01-01T00:00:00.000000000',
'2022-02-01T00:00:00.000000000', '2022-03-01T00:00:00.000000000'],
dtype='datetime64[ns]')[32344380 values with dtype=float32]
sst = file.sst.groupby(file.time.dt.year).mean().rename({'year':'time'})
sst = sst.sel(time=slice('1979','2021'))
sst
<xarray.DataArray 'sst' (time: 43, lat: 89, lon: 180)>
array([[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]]], dtype=float32)
Coordinates:
* lat (lat) float32 88.0 86.0 84.0 82.0 80.0 ... -82.0 -84.0 -86.0 -88.0
* lon (lon) float32 0.0 2.0 4.0 6.0 8.0 ... 350.0 352.0 354.0 356.0 358.0
* time (time) int64 1979 1980 1981 1982 1983 ... 2017 2018 2019 2020 2021array([[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]]], dtype=float32)array([ 88., 86., 84., 82., 80., 78., 76., 74., 72., 70., 68., 66.,
64., 62., 60., 58., 56., 54., 52., 50., 48., 46., 44., 42.,
40., 38., 36., 34., 32., 30., 28., 26., 24., 22., 20., 18.,
16., 14., 12., 10., 8., 6., 4., 2., 0., -2., -4., -6.,
-8., -10., -12., -14., -16., -18., -20., -22., -24., -26., -28., -30.,
-32., -34., -36., -38., -40., -42., -44., -46., -48., -50., -52., -54.,
-56., -58., -60., -62., -64., -66., -68., -70., -72., -74., -76., -78.,
-80., -82., -84., -86., -88.], dtype=float32)array([ 0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20., 22.,
24., 26., 28., 30., 32., 34., 36., 38., 40., 42., 44., 46.,
48., 50., 52., 54., 56., 58., 60., 62., 64., 66., 68., 70.,
72., 74., 76., 78., 80., 82., 84., 86., 88., 90., 92., 94.,
96., 98., 100., 102., 104., 106., 108., 110., 112., 114., 116., 118.,
120., 122., 124., 126., 128., 130., 132., 134., 136., 138., 140., 142.,
144., 146., 148., 150., 152., 154., 156., 158., 160., 162., 164., 166.,
168., 170., 172., 174., 176., 178., 180., 182., 184., 186., 188., 190.,
192., 194., 196., 198., 200., 202., 204., 206., 208., 210., 212., 214.,
216., 218., 220., 222., 224., 226., 228., 230., 232., 234., 236., 238.,
240., 242., 244., 246., 248., 250., 252., 254., 256., 258., 260., 262.,
264., 266., 268., 270., 272., 274., 276., 278., 280., 282., 284., 286.,
288., 290., 292., 294., 296., 298., 300., 302., 304., 306., 308., 310.,
312., 314., 316., 318., 320., 322., 324., 326., 328., 330., 332., 334.,
336., 338., 340., 342., 344., 346., 348., 350., 352., 354., 356., 358.],
dtype=float32)array([1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990,
1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014,
2015, 2016, 2017, 2018, 2019, 2020, 2021], dtype=int64)# eof
solver = Eof(sst)
eof = solver.eofs(neofs=3,eofscaling=2)
pc = solver.pcs(npcs=3,pcscaling=1)
print('eof:',eof.shape)
print('-------------------')
print('pc:',pc.shape)
eof: (3, 89, 180) ------------------- pc: (43, 3)
# 方差贡献率
var=solver.varianceFraction(neigs=3)
var
<xarray.DataArray 'variance_fractions' (mode: 3)>
array([0.32112256, 0.16426052, 0.06989542], dtype=float32)
Coordinates:
* mode (mode) int32 0 1 2
Attributes:
long_name: variance_fractionsarray([0.32112256, 0.16426052, 0.06989542], dtype=float32)
array([0, 1, 2])
year = np.arange(1979,2022)
def mapart(ax): #定义一个绘制地图地图的函数
'''
添加地图元素
'''
projection = ccrs.PlateCarree(central_longitude=180) #指定投影为经纬度投影,并指定中心经度为180°
ax.coastlines(color='k', lw=0.5) #添加海岸线
ax.add_feature(cfeature.LAND, facecolor='white') #添加陆地
# 设置地图范围,经度为(0,358),维度为(-88,88)
ax.set_extent([0, 358, -88, 88], crs=ccrs.PlateCarree(central_longitude=180))
# 设置经纬度标签
ax.set_xticks([-120, -60, 0, 60, 120], crs=projection)
ax.set_yticks([-80, -40, 0, 40, 80], crs=projection)
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
def eof_contourf(EOFs, PCs, pers, year):
plt.close
pers=(pers*100).values
fig = plt.figure(figsize=(18, 12), dpi=300)
projection = ccrs.PlateCarree(central_longitude=180)
ax1 = fig.add_subplot(3, 2, 1, projection=projection)
mapart(ax1)
p = ax1.contourf(file.lon,
file.lat,
EOFs[0],
np.linspace(-1,
1,
41),
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree())
ax1.set_title('mode1 (%s' % (round(pers[0], 2))+"%)", loc='left')
ax2 = fig.add_subplot(3, 2, 2)
ax2.plot(year, PCs[:, 0], color='k', linewidth=1.2, linestyle='--')
b = ax2.bar(year, PCs[:, 0], color='r')
for bar, height in zip(b, PCs[:, 0]):
if height < 0:
bar.set(color='blue')
ax2.set_title('PC1' % (round(pers[0], 2)), loc='left')
ax3 = fig.add_subplot(3, 2, 3, projection=projection)
mapart(ax3)
pp = ax3.contourf(file.lon,
file.lat,
EOFs[1],
np.linspace(-1,
1,
41),
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree())
ax3.set_title('mode2 (%s' % (round(pers[1], 2))+"%)", loc='left')
ax4 = fig.add_subplot(3, 2, 4)
ax4.plot(year, PCs[:, 1], color='k', linewidth=1.2, linestyle='--')
ax4.set_title('PC2' % (round(pers[1], 2)), loc='left')
bb = ax4.bar(year, PCs[:, 1], color='r')
for bar, height in zip(bb, PCs[:, 1]):
if height < 0:
bar.set(color='blue')
ax5 = fig.add_subplot(3, 2, 5, projection=projection)
mapart(ax5)
ppp = ax5.contourf(file.lon,
file.lat,
EOFs[2],
np.linspace(-1,
1,
41),
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree())
ax5.set_title('mode3 (%s' % (round(pers[2], 2))+"%)", loc='left')
ax6 = fig.add_subplot(3, 2, 6)
ax6.plot(year, PCs[:, 2], color='k', linewidth=1.2, linestyle='--')
ax6.set_title('PC3' % (round(pers[2], 2)), loc='left')
bbb = ax6.bar(year, PCs[:, 2], color='r')
for bar, height in zip(bbb, PCs[:, 2]):
if height < 0:
bar.set(color='blue')
# 为柱状图添加0标准线
ax2.axhline(y=0, linewidth=1, color='k', linestyle='-')
ax4.axhline(y=0, linewidth=1, color='k', linestyle='-')
ax6.axhline(y=0, linewidth=1, color='k', linestyle='-')
# 在图下边留白边放colorbar
fig.subplots_adjust(bottom=0.15)
# colorbar位置: 左 下 宽 高
l = 0.25
b = 0.04
w = 0.6
h = 0.015
# 对应 l,b,w,h;设置colorbar位置;
rect = [l, b, w, h]
cbar_ax = fig.add_axes(rect)
c = plt.colorbar(
pp,
cax=cbar_ax,
orientation='horizontal',
aspect=20,
pad=0.1) #绘制colorbar
c.ax.tick_params(labelsize=14) #设置colorbar的标签大小
plt.subplots_adjust(wspace=-0.1, hspace=0.3)
plt.show()
eof_contourf(eof*(-1),pc*(-1),var,year)
D:\anaconda\lib\site-packages\cartopy\crs.py:245: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(multi_line_string) > 1: D:\anaconda\lib\site-packages\cartopy\crs.py:297: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. for line in multi_line_string: D:\anaconda\lib\site-packages\cartopy\crs.py:364: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(p_mline) > 0: D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. line_strings = list(multi_line_string) D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. line_strings = list(multi_line_string) D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. line_strings.extend(multi_line_string) D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. line_strings.extend(multi_line_string)
north = solver.northTest(neigs=3).data
eig = solver.eigenvalues(neigs=3).data
print("north误差",north)
print("特征值",eig)
north误差 [115.06639 58.858727 25.045311] 特征值 [533.5409 272.9167 116.130325]
def northplt(north,eig,k):
ispass = False
for i in range (0,k-1):
if (eig[i+1]+north[i+1])<(eig[i]-north[i]):
print('第%d模态和第%d模态显著分离,通过north检验'%(i+1,i+2))
else:
print('第%d模态和第%d模态分离不显著,没有通过north检验'%(i+1,i+2))
break
else:
ispass = True
if ispass:
plt.errorbar(np.arange(k)+1, eig,yerr=north)
plt.xticks(np.arange(k)+1)
plt.title('north test PASS')
plt.show()
northplt(north,eig,3)
第1模态和第2模态显著分离,通过north检验 第2模态和第3模态显著分离,通过north检验
可以看到,三个主分量无重叠部分,显著分离
from scipy.stats import pearsonr
def nonan(a):
a = np.array(a)
mask = np.isnan(a)
a[mask] = 0
b = np.reshape(a,(89*180))
return b
def xiangguantest(eof,pc,k):
for i in range(k-1):
for j in range(i+1,k):
r1,p=pearsonr(nonan(eof[i,:,:]),nonan(eof[j,:,:]))
r1 = round(r1,2)
print(f'第{i+1}模态和第{j+1}模态空间相关系数:{r1}')
r2,p=pearsonr(pc[:,i],pc[:,i+1])
r2 = round(r2,2)
print(f'第{i+1}模态和第{j+1}模态时间相关系数:{r2}')
xiangguantest(eof,pc,3)
第1模态和第2模态空间相关系数:0.17 第1模态和第2模态时间相关系数:0.0 第1模态和第3模态空间相关系数:0.04 第1模态和第3模态时间相关系数:0.0 第2模态和第3模态空间相关系数:-0.02 第2模态和第3模态时间相关系数:-0.0
def mapart_(ax): #定义一个绘制地图地图的函数
'''
添加地图元素
'''
projection = ccrs.PlateCarree(central_longitude=180) #指定投影为经纬度投影,并指定中心经度为180°
ax.coastlines(color='k', lw=0.5) #添加海岸线
ax.add_feature(cfeature.LAND, facecolor='white') #添加陆地
# 设置地图范围,经度为(0,358),维度为(-88,88)
ax.set_extent([0, 358, -88, 88], crs=ccrs.PlateCarree(central_longitude=180))
# 设置经纬度标签
ax.set_xticks([-120, -60, 0, 60, 120], crs=projection)
ax.set_yticks([-80, -40, 0, 40, 80], crs=projection)
lon_formatter = LongitudeFormatter(zero_direction_label=True) #给标签添加对应的N,S,E,W
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
def eof_contourf_(EOFs, PCs, pers, year): #定义一个绘制填充图的函数
plt.close
pers=(pers*100).values #将方差转换为百分数的形式
fig = plt.figure(figsize=(18, 12), dpi=300)
projection = ccrs.PlateCarree(central_longitude=180)
ax1 = fig.add_subplot(3, 2, 1, projection=projection)
mapart(ax1)
p = ax1.contourf(file.lon,
file.lat,
EOFs[0],
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree())
ax1.set_title('mode1 (%s' % (round(pers[0], 2))+"%)", loc='left')
ax2 = fig.add_subplot(3, 2, 2) #绘制第二个子图
ax2.plot(year, PCs[:, 0], color='k', linewidth=1.2, linestyle='--')
b = ax2.bar(year, PCs[:, 0], color='r')
for bar, height in zip(b, PCs[:, 0]):
if height < 0:
bar.set(color='blue')
ax2.set_title('PC1' % (round(pers[0], 2)), loc='left')
ax3 = fig.add_subplot(3, 2, 3, projection=projection)
mapart(ax3)
pp = ax3.contourf(file.lon,
file.lat,
EOFs[1],
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree())
ax3.set_title('mode2 (%s' % (round(pers[1], 2))+"%)", loc='left')
ax4 = fig.add_subplot(3, 2, 4)
ax4.plot(year, PCs[:, 1], color='k', linewidth=1.2, linestyle='--')
ax4.set_title('PC2' % (round(pers[1], 2)), loc='left')
bb = ax4.bar(year, PCs[:, 1], color='r')
for bar, height in zip(bb, PCs[:, 1]):
if height < 0:
bar.set(color='blue')
ax5 = fig.add_subplot(3, 2, 5, projection=projection)
mapart(ax5)
ppp = ax5.contourf(file.lon,
file.lat,
EOFs[2],
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree())
ax5.set_title('mode3 (%s' % (round(pers[2], 2))+"%)", loc='left')
ax6 = fig.add_subplot(3, 2, 6)
ax6.plot(year, PCs[:, 2], color='k', linewidth=1.2, linestyle='--')
ax6.set_title('PC3' % (round(pers[2], 2)), loc='left')
bbb = ax6.bar(year, PCs[:, 2], color='r')
for bar, height in zip(bbb, PCs[:, 2]):
if height < 0:
bar.set(color='blue')
# 为柱状图添加0标准线
ax2.axhline(y=0, linewidth=1, color='k', linestyle='-')
ax4.axhline(y=0, linewidth=1, color='k', linestyle='-')
ax6.axhline(y=0, linewidth=1, color='k', linestyle='-')
# 在图下边留白边放colorbar
fig.subplots_adjust(bottom=0.15)
# colorbar位置: 左 下 宽 高
l = 0.25
b = 0.04
w = 0.6
h = 0.015
# 对应 l,b,w,h;设置colorbar位置;
rect = [l, b, w, h]
cbar_ax = fig.add_axes(rect)
c = plt.colorbar(
pp,
cax=cbar_ax,
orientation='horizontal',
aspect=20,
pad=0.1) #绘制colorbar
c.ax.tick_params(labelsize=14) #设置colorbar的标签大小
plt.subplots_adjust(wspace=-0.1, hspace=0.3)
# plt.tight_layout()
# plt.savefig('./figure/eof2.png')
plt.show()
# eof
solver = Eof(sst,center=False)
# solver = Eof(sst)
eof_ = solver.eofs(neofs=3,eofscaling=2)
pc_ = solver.pcs(npcs=3,pcscaling=1)
var_=solver.varianceFraction(neigs=3)
print('eof_:',eof_.shape)
print('-------------------')
print('pc_:',pc_.shape)
eof_: (3, 89, 180) ------------------- pc_: (43, 3)
eof_contourf_(eof_*(-1),pc_*(-1),var_,year)
north_ = solver.northTest(neigs=3).data
eig_ = solver.eigenvalues(neigs=3).data
print("north误差",north_)
print("特征值",eig_)
north误差 [7.8232719e+05 8.7714035e+01 4.2884262e+01] 特征值 [3.6275018e+06 4.0671323e+02 1.9884615e+02]
northplt(north_,eig_,3)
第1模态和第2模态显著分离,通过north检验 第2模态和第3模态显著分离,通过north检验
xiangguantest(eof_,pc_,3)
第1模态和第2模态空间相关系数:0.1 第1模态和第2模态时间相关系数:-0.59 第1模态和第3模态空间相关系数:0.12 第1模态和第3模态时间相关系数:-0.59 第2模态和第3模态空间相关系数:-0.02 第2模态和第3模态时间相关系数:-0.0
from sklearn import preprocessing
import eofs
sst_scaled_t = np.zeros((43,89,180))
for i in range(43):
sst_scaled_t[i,:,:] = preprocessing.scale(sst.data[i,:,:])
sst_scaled_t
D:\anaconda\lib\site-packages\sklearn\preprocessing\_data.py:235: UserWarning: Numerical issues were encountered when centering the data and might not be solved. Dataset may contain too large values. You may need to prescale your features. warnings.warn( D:\anaconda\lib\site-packages\sklearn\preprocessing\_data.py:254: UserWarning: Numerical issues were encountered when scaling the data and might not be solved. The standard deviation of the data is probably very close to 0. warnings.warn(
array([[[-1.32081151, -1.31283498, -1.35773385, ..., -1.27647579,
-1.29350352, -1.29724193],
[-1.32081151, -1.31283498, -1.35773385, ..., -1.27647579,
-1.29350352, -1.29724193],
[-1.32081151, -1.31283498, -1.35773385, ..., -1.27647579,
-1.29350352, -1.29724193],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.33402646, -1.32488191, -1.36965883, ..., -1.29251158,
-1.30991757, -1.31138611],
[-1.33402646, -1.32488191, -1.36965883, ..., -1.29251158,
-1.30991757, -1.31138611],
[-1.33402646, -1.32488191, -1.36965883, ..., -1.29251158,
-1.30991757, -1.31138611],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.32527184, -1.31729436, -1.36138654, ..., -1.28340483,
-1.29943419, -1.30155838],
[-1.32527184, -1.31729436, -1.36138654, ..., -1.28340483,
-1.29943419, -1.30155838],
[-1.32527184, -1.31729436, -1.36138654, ..., -1.28340483,
-1.29943419, -1.30155838],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
...,
[[-1.35812294, -1.35505021, -1.40585339, ..., -1.31424212,
-1.33291507, -1.33418429],
[-1.35812294, -1.35505021, -1.40585339, ..., -1.31424212,
-1.33291507, -1.33418429],
[-1.35812294, -1.35505021, -1.40585339, ..., -1.31424212,
-1.33291507, -1.33418429],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.35244596, -1.35033238, -1.39882708, ..., -1.30525815,
-1.32402062, -1.32561731],
[-1.35244596, -1.35033238, -1.39882708, ..., -1.30525815,
-1.32402062, -1.32561731],
[-1.35244596, -1.35033238, -1.39882708, ..., -1.30525815,
-1.32402062, -1.32561731],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.35440552, -1.34971094, -1.39754999, ..., -1.31392741,
-1.3279072 , -1.32902801],
[-1.35440552, -1.34971094, -1.39754999, ..., -1.31392741,
-1.3279072 , -1.32902801],
[-1.35440552, -1.34971094, -1.39754999, ..., -1.31392741,
-1.3279072 , -1.32902801],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]]])
solver = eofs.standard.Eof(sst_scaled_t)
eof__ = solver.eofs(neofs=3,eofscaling=2)
pc__ = solver.pcs(npcs=3,pcscaling=1)
var__=solver.varianceFraction(neigs=3)
print('eof__:',eof__.shape)
print('-------------------')
print('pc__:',pc__.shape)
eof__: (3, 89, 180) ------------------- pc__: (43, 3)
eof_contourf_(eof__*(-1),pc__*(-1),var_,year)
D:\anaconda\lib\site-packages\cartopy\crs.py:245: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(multi_line_string) > 1: D:\anaconda\lib\site-packages\cartopy\crs.py:297: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. for line in multi_line_string: D:\anaconda\lib\site-packages\cartopy\crs.py:364: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(p_mline) > 0: D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. line_strings = list(multi_line_string) D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. line_strings = list(multi_line_string) D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. line_strings.extend(multi_line_string) D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. line_strings.extend(multi_line_string)
north__ = solver.northTest(neigs=3)
eig__ = solver.eigenvalues(neigs=3)
print("north误差",north__)
print("特征值",eig__)
north误差 [0.51234382 0.24994451 0.15094646] 特征值 [2.37564056 1.15894504 0.69990995]
northplt(north__,eig__,3)
第1模态和第2模态显著分离,通过north检验 第2模态和第3模态显著分离,通过north检验
xiangguantest(eof__,pc__,3)
第1模态和第2模态空间相关系数:-0.0 第1模态和第2模态时间相关系数:0.0 第1模态和第3模态空间相关系数:0.0 第1模态和第3模态时间相关系数:0.0 第2模态和第3模态空间相关系数:0.0 第2模态和第3模态时间相关系数:0.0
sst
<xarray.DataArray 'sst' (time: 43, lat: 89, lon: 180)>
array([[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]]], dtype=float32)
Coordinates:
* lat (lat) float32 88.0 86.0 84.0 82.0 80.0 ... -82.0 -84.0 -86.0 -88.0
* lon (lon) float32 0.0 2.0 4.0 6.0 8.0 ... 350.0 352.0 354.0 356.0 358.0
* time (time) int64 1979 1980 1981 1982 1983 ... 2017 2018 2019 2020 2021array([[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]],
[[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
[-1.7999998, -1.7999998, -1.7999998, ..., -1.7999998,
-1.7999998, -1.7999998],
...,
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan],
[ nan, nan, nan, ..., nan,
nan, nan]]], dtype=float32)array([ 88., 86., 84., 82., 80., 78., 76., 74., 72., 70., 68., 66.,
64., 62., 60., 58., 56., 54., 52., 50., 48., 46., 44., 42.,
40., 38., 36., 34., 32., 30., 28., 26., 24., 22., 20., 18.,
16., 14., 12., 10., 8., 6., 4., 2., 0., -2., -4., -6.,
-8., -10., -12., -14., -16., -18., -20., -22., -24., -26., -28., -30.,
-32., -34., -36., -38., -40., -42., -44., -46., -48., -50., -52., -54.,
-56., -58., -60., -62., -64., -66., -68., -70., -72., -74., -76., -78.,
-80., -82., -84., -86., -88.], dtype=float32)array([ 0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20., 22.,
24., 26., 28., 30., 32., 34., 36., 38., 40., 42., 44., 46.,
48., 50., 52., 54., 56., 58., 60., 62., 64., 66., 68., 70.,
72., 74., 76., 78., 80., 82., 84., 86., 88., 90., 92., 94.,
96., 98., 100., 102., 104., 106., 108., 110., 112., 114., 116., 118.,
120., 122., 124., 126., 128., 130., 132., 134., 136., 138., 140., 142.,
144., 146., 148., 150., 152., 154., 156., 158., 160., 162., 164., 166.,
168., 170., 172., 174., 176., 178., 180., 182., 184., 186., 188., 190.,
192., 194., 196., 198., 200., 202., 204., 206., 208., 210., 212., 214.,
216., 218., 220., 222., 224., 226., 228., 230., 232., 234., 236., 238.,
240., 242., 244., 246., 248., 250., 252., 254., 256., 258., 260., 262.,
264., 266., 268., 270., 272., 274., 276., 278., 280., 282., 284., 286.,
288., 290., 292., 294., 296., 298., 300., 302., 304., 306., 308., 310.,
312., 314., 316., 318., 320., 322., 324., 326., 328., 330., 332., 334.,
336., 338., 340., 342., 344., 346., 348., 350., 352., 354., 356., 358.],
dtype=float32)array([1979, 1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990,
1991, 1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002,
2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013, 2014,
2015, 2016, 2017, 2018, 2019, 2020, 2021], dtype=int64)
sst_scaled_p = np.zeros((43,89,180))
for i in range(89):
for j in range(180):
sst_scaled_p[:,i,j] = preprocessing.scale(sst.data[:,i,j])
sst_scaled_p
D:\anaconda\lib\site-packages\sklearn\preprocessing\_data.py:235: UserWarning: Numerical issues were encountered when centering the data and might not be solved. Dataset may contain too large values. You may need to prescale your features. warnings.warn( IOPub data rate exceeded. The Jupyter server will temporarily stop sending output to the client in order to avoid crashing it. To change this limit, set the config variable `--ServerApp.iopub_data_rate_limit`. Current values: ServerApp.iopub_data_rate_limit=1000000.0 (bytes/sec) ServerApp.rate_limit_window=3.0 (secs)
solver = eofs.standard.Eof(sst_scaled_p)
eof__p = solver.eofs(neofs=3,eofscaling=2)
pc__p = solver.pcs(npcs=3,pcscaling=1)
var__p=solver.varianceFraction(neigs=3)
print('eof__p:',eof__p.shape)
print('-------------------')
print('pc__p:',pc__p.shape)
eof__p: (3, 89, 180) ------------------- pc__p: (43, 3)
def eof_contourf__(EOFs, PCs, pers, year): #定义一个绘制填充图的函数
plt.close
pers=(pers*100) #将方差转换为百分数的形式
fig = plt.figure(figsize=(18, 12), dpi=300)
projection = ccrs.PlateCarree(central_longitude=180)
ax1 = fig.add_subplot(3, 2, 1, projection=projection)
mapart(ax1)
p = ax1.contourf(file.lon,
file.lat,
EOFs[0],
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree())
ax1.set_title('mode1 (%s' % (round(pers[0], 2))+"%)", loc='left')
ax2 = fig.add_subplot(3, 2, 2) #绘制第二个子图
ax2.plot(year, PCs[:, 0], color='k', linewidth=1.2, linestyle='--')
b = ax2.bar(year, PCs[:, 0], color='r')
for bar, height in zip(b, PCs[:, 0]):
if height < 0:
bar.set(color='blue')
ax2.set_title('PC1' % (round(pers[0], 2)), loc='left')
ax3 = fig.add_subplot(3, 2, 3, projection=projection)
mapart(ax3)
pp = ax3.contourf(file.lon,
file.lat,
EOFs[1],
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree())
ax3.set_title('mode2 (%s' % (round(pers[1], 2))+"%)", loc='left')
ax4 = fig.add_subplot(3, 2, 4)
ax4.plot(year, PCs[:, 1], color='k', linewidth=1.2, linestyle='--')
ax4.set_title('PC2' % (round(pers[1], 2)), loc='left')
bb = ax4.bar(year, PCs[:, 1], color='r')
for bar, height in zip(bb, PCs[:, 1]):
if height < 0:
bar.set(color='blue')
ax5 = fig.add_subplot(3, 2, 5, projection=projection)
mapart(ax5)
ppp = ax5.contourf(file.lon,
file.lat,
EOFs[2],
cmap=cmaps.BlueWhiteOrangeRed,
transform=ccrs.PlateCarree())
ax5.set_title('mode3 (%s' % (round(pers[2], 2))+"%)", loc='left')
ax6 = fig.add_subplot(3, 2, 6)
ax6.plot(year, PCs[:, 2], color='k', linewidth=1.2, linestyle='--')
ax6.set_title('PC3' % (round(pers[2], 2)), loc='left')
bbb = ax6.bar(year, PCs[:, 2], color='r')
for bar, height in zip(bbb, PCs[:, 2]):
if height < 0:
bar.set(color='blue')
# 为柱状图添加0标准线
ax2.axhline(y=0, linewidth=1, color='k', linestyle='-')
ax4.axhline(y=0, linewidth=1, color='k', linestyle='-')
ax6.axhline(y=0, linewidth=1, color='k', linestyle='-')
# 在图下边留白边放colorbar
fig.subplots_adjust(bottom=0.15)
# colorbar位置: 左 下 宽 高
l = 0.25
b = 0.04
w = 0.6
h = 0.015
# 对应 l,b,w,h;设置colorbar位置;
rect = [l, b, w, h]
cbar_ax = fig.add_axes(rect)
c = plt.colorbar(
pp,
cax=cbar_ax,
orientation='horizontal',
aspect=20,
pad=0.1) #绘制colorbar
c.ax.tick_params(labelsize=14) #设置colorbar的标签大小
plt.subplots_adjust(wspace=-0.1, hspace=0.3)
# plt.tight_layout()
# plt.savefig('./figure/eof2.png')
plt.show()
eof_contourf__(eof__p*(-1),pc__p*(-1),var__p,year)
D:\anaconda\lib\site-packages\cartopy\crs.py:245: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(multi_line_string) > 1: D:\anaconda\lib\site-packages\cartopy\crs.py:297: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. for line in multi_line_string: D:\anaconda\lib\site-packages\cartopy\crs.py:364: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(p_mline) > 0: D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. line_strings = list(multi_line_string) D:\anaconda\lib\site-packages\cartopy\crs.py:256: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. line_strings = list(multi_line_string) D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. line_strings.extend(multi_line_string) D:\anaconda\lib\site-packages\cartopy\crs.py:402: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. line_strings.extend(multi_line_string)
north__p = solver.northTest(neigs=3)
eig__p = solver.eigenvalues(neigs=3)
print("north误差",north__p)
print("特征值",eig__p)
north误差 [686.41866511 301.32724554 152.68502965] 特征值 [3182.79241419 1397.1969587 707.97135746]
northplt(north__,eig__,3)
第1模态和第2模态显著分离,通过north检验 第2模态和第3模态显著分离,通过north检验
xiangguantest(eof__p,pc__p,3)
第1模态和第2模态空间相关系数:0.2 第1模态和第2模态时间相关系数:0.0 第1模态和第3模态空间相关系数:-0.01 第1模态和第3模态时间相关系数:0.0 第2模态和第3模态空间相关系数:0.01 第2模态和第3模态时间相关系数:-0.0